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1.  Introduction 

Despite  significant  progress  in  computational  sciences,  challenges  persist  in 
the  accurate  numerical  simulation  of  a broad  spectrum  of  dynamic,  multi- 
physics phenomena  relevant  to  aerospace  systems.  These  challenging  areas 
include  the  direct  and  large-eddy  simulation  of  turbulence,  aeroacoustics, 
fluid/structure  interactions,  electromagnetics,  and  magneto-gasdynamics. 
In  order  to  reduce  the  severe  computational  requirements  of  standard  low- 
order  schemes,  higher-order  formulations,  as  well  as  massively  parallel  ap- 
proaches are  being  actively  pursued.  Due  to  their  spectral-like  resolution 
and  ease  of  extension  to  multiple  disciplines,  high-order  compact  schemes[l] 
represent  an  attractive  choice  for  reducing  dispersion,  anisotropy  and  dis- 
sipation errors  associated  with  spatial  discretizations.  Until  recently,  these 
schemes  have  mostly  been  used  in  conjunction  with  explicit  time-integration 
methods  to  address  complex  flow  physics  on  simple  Cartesian  geometries. 

Recent  work  performed  at  the  Air  Force  Research  Laboratory  [2,  3,  4, 
5,  6,  7]  has  extended  the  use  of  compact  algorithms  to  more  practical  ap- 
plications. This  has  been  achieved  through  the  development  and  improved 
treatment  of  the  various  critical  elements  comprising  the  overall  numerical 
approach.  Particular  attention  has  been  focused  on  enhanced  high-order 
(up  to  lOth-order)  low-pass  filtering  techniques,  accurate  and  robust  near- 
boundary formulations,  proper  metric  evaluation,  multi-domain  implemen- 
tation strategies,  and  sub-iterative  implicit  time-advancement  methods.  As 
an  outcome  of  this  sustained  effort,  the  highly  accurate  compact  algorithm 
has  been  successfully  applied  to  efficiently  solve  a range  of  multi-physics 
phenomena  described  by  the  Euler,  Navier-Stokes  and  MHD  equations  on 
3-D  curvilinear  and  dynamic  grids  using  either  explicit  or  implicit  time 
integration  approaches. 
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These  high-fidelity  computational  tools  are  currently  being  transitioned 
to  the  multidisciplinary  simulation  of  complex  phenomena  relevant  to  Air 
Force  systems,  including:  weapon-bay  cavity  acoustics,  hypersonic  flow  con- 
trol, high-  angle-of-attack  aerodynamics,  and  non-linear  aeroelastic  response. 
A brief  description  of  the  governing  equations,  the  various  elements  of  the 
numerical  approach,  as  well  as  a few  representative  applications  are  in- 
cluded in  this  paper. 


2.  Governing  Equations 


In  order  to  develop  a procedure  suitable  for  nonlinear  fluid  dynamic,  aeroa- 
coustic  and  aeroelastic  applications  over  complex  configurations,  the  full 
Navier-Stokes  equations  are  selected  and  are  cast  in  strong  conservative 
form  introducing  a general  time-dependent  curvilinear  coordinate  trans- 
formation ( x,y,z,t ) — » (£,  rj,  (,  r).  In  vector  notation,  and  employing  non- 
dimensional  variables,  these  equations  are: 


d_ 

dr 


OF  dG  dH 
+ d£  + dr,  + dC  ~ 


1 9FV  3GV 

Re[  d£  dr. 


dHv 

dC 


} + S/J 


(1) 


where  U = {p,  pu , pv,  pw,  pE}  denotes  the  solution  vector,  F,  G,  H,  FV,GV, 
Hv  are  the  fluxes,  S denotes  a source  term,  and  J = d (£,  rj,  £,  r)  /d  (x,  y,  z,  t) 
is  the  transformation  Jacobian  which  for  dynamic  meshes  is  a function  of 
time. 


3.  Numerical  Method 
3.1.  SPATIAL  DISCRETIZATION 

A finite-difference  approach  is  employed  to  discretize  the  governing  equa- 
tions. This  choice  is  motivated  by  the  relative  ease  of  formal  extension  to 
higher-order  accuracy.  For  any  scalar  quantity,  <p,  such  as  a metric,  flux 
component  or  flow  variable,  the  spatial  derivative  $ is  obtained  in  the 
transformed  plane  by  solving  the  tridiagonal  system: 

Pi/  i i/  i pi/  tpi—2  | 1 1 /o', 

r^i-i  + + r<£i+1  = 6 +a (2) 

where  T,  a and  b determine  the  spatial  properties  of  the  algorithm.  The  for- 
mula yields  the  compact  five-point,  sixth-order  C6 , and  three-point  fourth- 
order  C4  schemes  [1]  with  r = 1/3,  a = 14/9,  b — 1/9  and  T = 1/4, 
a = 3/2,  6 = 0 respectively.  Equation  (2)  also  incorporates  the  standard 
explicit  fourth-order  ( E4 ) and  second-order  (E2)  approaches  for  which  the 
coefficients  are  (T  = 0,  a = 4/3,  b = —1/3)  and  (T  = 0,  a = 1,  6 = 0)  re- 
spectively. The  dispersion  characteristics  and  truncation  error  of  the  above 
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schemes  can  be  found  in  Refs.  [1,  5].  It  should  be  noted  that  for  a given 
order  of  accuracy,  the  compact  schemes  are  significantly  superior  to  their 
explicit  (non-compact)  counterparts.  The  scheme  of  Eqn.  (2)  cannot  be  ap- 
plied at  a few  points  near  each  boundary  where  the  stencil  protrudes  the 
domain.  Here,  special  one-sided  boundary  schemes,  such  as  those  described 
in  Refs.  [5,  8]  are  employed. 

In  order  to  compute  the  residual,  the  derivatives  of  the  inviscid  fluxes  are 
obtained  by  first  forming  the  fluxes  at  the  nodes  and  subsequently  differenti- 
ating each  component  with  the  above  formulas.  Viscous  terms  are  obtained 
by  first  computing  derivatives  of  the  primitive  variables.  The  components 
of  the  viscous  flux  are  then  constructed  at  each  node  and  differentiated  by 
a second  application  of  the  same  scheme.  Although  this  approach  is  not 
as  accurate  as  that  in  which  a Pade-type  scheme  is  employed  directly  for 
the  second-derivative,  it  is  significantly  cheaper  to  implement  in  curvilinear 
coordinates.  As  previously  demonstrated  in  Ref.  [2],  successive  differentia- 
tion yields  an  accurate  and  stable  method  in  conjunction  with  the  added 
filtering  component  which  is  described  next. 

3.2.  HIGH-ORDER  FILTERING  SCHEME 

Compact-difference  discretizations,  like  other  centered  schemes,  are  non- 
dissipative  and  are  therefore  susceptible  to  numerical  instabilities  due  to  the 
unrestricted  growth  of  high-frequency  modes.  These  difficulties  originate 
from  several  sources  including  mesh  non-uniformity,  approximate  bound- 
ary conditions  and  nonlinear  flow  features.  In  order  to  extend  the  present 
solver  to  practical  simulations,  while  retaining  the  improved  accuracy  of  the 
spatial  compact  discretization,  a high-order  implicit  filtering  technique  [2,  4] 
is  incorporated.  If  a typical  component  of  the  solution  vector  is  denoted  by 
(j>,  filtered  values  <f>  satisfy, 

1 T 4>i  T CXf<j>i+ 1 = Sra_Q-^-  ( (f>i+n  + 4>i—n ) (3) 

Equation  (3),  with  the  proper  choice  of  coefficients,  provides  an  21Vth-order 
formula  on  a 2N  + 1 point  stencil.  The  N + 1 coefficients,  aQ,  oi, a^, 
are  derived  in  terms  of  ay  with  Taylor-  and  Fourier-series  analyses  and 
are  found  in  Refs.  [2,  3]  along  with  some  spectral  filter  responses.  The  ad- 
justable parameter  ay  satisfies  the  inequality  —0.5  < ay  < 0.5,  with  higher 
values  of  ay  corresponding  to  a less  dissipative  filter.  In  multi-dimensional 
problems  the  filter  operator  is  applied  sequentially  in  each  coordinate  di- 
rection. For  the  near-boundary  points,  the  filtering  strategies  described  in 
Refs.  [2,  3]  are  employed.  Up  to  tenth-order  filter  formulas  have  been  suc- 
cessfully applied  to  solve  the  Maxwell  [4],  Navier-Stokes  [2,  9,  10,  11],  and 
MHD  [12]  equations  in  curvilinear  geometries. 


16  M.  VISBAL,  D.  GAITONDE  AND  D.  RIZZETTA 

3.3.  METRIC  EVALUATION  FOR  CURVILINEAR  DYNAMIC  MESHES 

The  extension  of  high-order  schemes  to  3-D  curvilinear  meshes  demands 
that  issues  of  freestream  preservation  and  metric  cancellation  be  carefully 
addressed.  These  errors,  which  arise  in  finite-difference  discretizations  of 
governing  equations  written  in  strong-conservation  form,  can  catastrophi- 
cally degrade  the  fidelity  of  standard  second-order  as  well  as  higher-order 
approaches  [3].  In  deriving  the  flow  equations  in  strong-conservation  form, 
the  following  metric  identities  have  been  implicitly  invoked, 


h = (Sx/JU  + (Vx/J)v  + (Cx/J)  c = 0 (4) 

h = (Zy/J)t  + (%/J),  + (C y/Jk  = 0 (5) 

h = (UJ)t  + (%/J)n  + (Cz/J)  C = 0 (6) 

h = (l/J)r  + + ( TH/J)v  + (Ct/J)  c = 0 (7) 


In  Ref.  [2]  it  was  shown  that  on  stretched  curvilinear  2-D  meshes,  the  com- 
pact scheme  exhibits  freestream  preservation  when  the  metrics  are  evalu- 
ated with  the  same  finite-difference  expressions  as  those  employed  for  the 
fluxes.  It  was  also  demonstrated  that  the  practice  of  prescribing  analytic 
metrics  on  stretched  curvilinear  grids  can  lead  to  unacceptable  errors  and 
therefore  should  in  general  be  avoided.  The  previous  straightforward  ap- 
proach of  calculating  the  metrics,  although  effective  in  2-D,  fails  to  provide 
metric  cancellation  for  general  3-D  curvilinear  configurations.  To  illustrate 
this  point,  consider  the  standard  metric  relations: 

^x/ J — Vrt^C,  — yc,zri 

Vx/J  = - V(,z<,  (8) 

C x/J  = Vizn  - ynz(, 

associated  with  the  identity  I\.  Evaluation  of  the  y and  z derivatives  in  the 
previous  expressions  using  explicit  or  compact  centered  schemes  does  not 
satisfy  I\,  and  as  a result,  significant  grid-induced  errors  may  appear.  To 
extend  the  high-order  compact  scheme  to  general  geometries,  the  metric 
terms  are  rewritten  prior  to  discretization  in  the  equivalent  (‘conservative’) 
form  [13]: 

ix/J  = (y»,*)c  - (y<;z)r) 

Vx/J  = (y<*)e  - (y^)<  (9) 

C x/J  = {y^z)t]  — (yT)z)z 

Similar  expressions  are  employed  for  the  remaining  metric  terms  in  order  to 
satisfy  the  identities  1 2 and  I3  above.  When  the  transformation  metrics  are 
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recast  in  this  manner,  and  the  derivatives  are  evaluated  with  the  same  high- 
order  formulas  employed  for  the  fluxes,  freestream  preservation  is  again 
recovered  in  general  time-invariant  3-D  curvilinear  geometries  [3]. 

For  deforming  and  moving  meshes,  the  identity  I4  must  be  also  satisfied 
to  eliminate  metric  cancellation  errors  and  to  ensure  freestream  preser- 
vation [6].  This  metric  identity  is  referred  to  in  the  literature  [13]  as  the 
Geometric  Conservation  Law  ( GCL ).  For  the  time-integration  methods  em- 
ployed in  this  work,  the  time-derivative  term  in  Eqn.  (1)  is  split  using 
chain-rule  differentiation  as  follows: 

(U/J)  r = {1/J)UT  + U(1/J)T  (10) 

Rather  than  attempting  to  compute  the  time  derivative  of  the  inverse  Ja- 
cobian directly  from  the  grid  coordinates  at  various  time  levels  (either  ana- 
lytically or  numerically),  we  simply  invoke  the  GCL  identity  I4  to  evaluate 


(1/J)T,  i.e. 

(1  /J)r 

— - [(£t/*7)f  + ivt/ J)t)  + (Ct/J)c] 

(11) 

where 

II 

-M^/J)+Vr^V/J)+Zr(UJ)] 

II 

ST 

~ [XT  iVx/J)  + VriVy/J)  + zt{Vz/ J)] 

(12) 

II 

-[Xt{(x/J)  +Vt{CV/J)  +Zt(Cz/J)\ 

For  the  case  of  an  analytically  prescribed  dynamic  mesh  transformation,  the 
grid  speeds  (xT,  yT , zT ) are  obtained  from  the  corresponding  analytic  expres- 
sions. An  example  in  which  the  grid  speeds  are  known  analytically  corre- 
sponds to  the  case  of  a maneuvering  wing  when  the  entire  mesh  moves  in 
a rigid  fashion.  In  many  practical  applications  involving  deforming  meshes 
(e.g.  dynamic  aeroelastic  simulations),  the  grid  speeds  are  not  known  a 
priori,  and  must  therefore  be  approximated  to  the  desired  degree  of  ac- 
curacy employing  the  evolving  grid  coordinates  at  several  time  levels.  As 
demonstrated  in  Ref.  [6],  the  high-order  method  retains  its  superior  accu- 
racy on  rapidly  distorting  meshes  when  the  procedure  outlined  above  is 
incorporated  for  the  time  metrics. 

3.4.  TIME-INTEGRATION  SCHEME 

Two  different  time-integration  approaches  are  incorporated  in  the  present 
family  of  solvers.  For  wave  propagation  applications,  the  equations  are  in- 
tegrated in  time  with  the  classical  fourth-order  four-stage  Runge-Kutta 
method  ( RK4 )•  The  scheme  is  implemented  in  low  storage  form  requiring 
three  levels  of  storage.  For  wall-bounded  viscous  flows,  the  stability  con- 
straint of  the  explicit  time-marching  scheme  is  found  to  render  the  approach 
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too  restrictive  and  inefficient.  Therefore,  the  implicit,  approximately-factored 
method  of  Beam  and  Warming  [14]  is  also  incorporated  and  augmented 
through  the  use  of  Newton-like  subiterations  in  order  to  achieve  second- 
order  time  accuracy.  In  delta  form,  the  scheme  may  be  written  as 


J 


_ip+i 


c(2)  ( dFP' 


JP+1X 


+ ft  At 8 £ , du 

J~1P+1  + ft  At8^  (#)l  Jp+1x 
J~1P+1+ft 
= —ft At  [ J-1P+1  (l+<t>)UP-'^+Wun+<t>un 

+UP{1/J)TP  + (fo)  + 

k ft)  + k (a')] 


(13) 


where  dF/dU  etc  are  flux  Jacobians,  8 represents  the  spatial  difference 
operator  and  AU  = Up+1  — Up.  For  improved  efficieny,  the  method  incor- 
porates the  diagonalization  procedure  of  Ref.  [15].  In  addition,  nonlinear 
artificial  dissipation  terms  [16]  are  appended  to  the  implicit  operator  to  en- 
hance stability.  Note  that  while  the  derivatives  of  the  flux  Jacobians  have 
been  obtained  to  second-order  accuracy  (denoted  with  the  superscript  (2)  in 
Eqn.  (13)),  those  on  the  right  hand  side,  i.e.,  in  the  residual,  are  evaluated 
with  the  compact-difference  higher-order  method.  In  order  to  reduce  errors 
associated  with  these  simplifications,  a sub-iteration  strategy  is  employed. 
Thus,  for  the  first  subiteration,  p = 1,  Up  = Un  and  as  p -»  oo,  Up  — > Un+1 . 
Typically,  three  subiterations  are  applied  per  time  step.  A range  of  numer- 
ical experiments  suggests  that  second-order  accuracy  in  time  is  adequate 
for  the  problems  considered  [6]. 


3.5.  MULTI-DOMAIN  STRATEGY 

Domain-decomposition  techniques  constitute  an  important  component  of 
modern  computational  strategy.  Due  to  their  spatially  implicit  nature, 
Pade-type  schemes  are  more  difficult  to  utilize  in  a multi-domain  environ- 
ment than  explicit  methods.  However,  a finite-size  overlap  can  be  employed 
with  the  present  compact /filtering  methodology  to  generate  a powerful  ap- 
proach applicable  to  complicated  curvilinear  meshes  [3].  Figure  1 depicts 
schematically  the  problem  of  a vortex  traveling  to  the  right  in  a recti- 
linear path.  The  domain  of  computation  is  divided  into  two  parts  to  be 
distributed  to  different  processors.  Each  sub-domain  is  supplemented  with 
several  points  from  the  adjacent  sub-domain  to  form  an  overlap  region, 
whose  details  for  a five-point  vertical  overlap  are  also  shown  in  Fig.  1. 
Although  the  overlap  points  are  collocated  they  have  been  shown  slightly 
staggered  for  clarity.  Each  vertical  line  is  denoted  by  its  i-index.  Data  is 
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exchanged  between  adjacent  subdomains  at  the  end  of  each  sub-iteration  of 
the  implicit  scheme  (or  each  stage  of  RK4 ),  as  well  as  after  each  application 
of  the  filter.  The  values  at  points  1 and  2 of  Mesh  2 are  set  to  be  identically 
equal  to  the  corresponding  updated  values  at  points  IL  — 4 and  IL  — 3 
of  Mesh  1.  Similarly,  reciprocal  information  is  transferred  through  points  4 
and  5 of  Mesh  2 which  “donate”  values  to  points  IL  — 1 and  IL  of  Mesh  1. 
More  details  on  the  accuracy  and  robustness  of  the  present  multi-domain 
approach  can  be  found  in  Ref.  [3]. 

4.  Results 

The  previous  computational  methodology  has  been  successfully  demon- 
strated for  a number  of  applications,  including:  unsteady  vortical  flows  [2], 
DNS  of  pulsed  walljet  [9]  and  synthetic  jet  actuators  [17],  LES  of  subsonic 
and  supersonic  flows  [10,  11],  non-linear  fluid  structure  interaction  [18], 
and  benchmark  problems  for  acoustic  scattering  [7]  and  MHD  [12].  In  the 
interest  of  brevity,  only  a subset  of  these  results  is  presented  here. 

4.1.  MULTI-DOMAIN  ACOUSTIC  SCATTERING  SIMULATION 

The  low  dispersion  error  characteristic  of  compact-difference  schemes  is  an 
attractive  property  in  the  simulation  of  wave  propagation  phenomena  asso- 
ciated with  acoustic  and  electromagnetic  scattering.  In  order  to  demon- 
strate the  capability  of  the  present  numerical  approach  to  treat  acous- 
tic phenomena  in  a multiple-domain  situation,  consider  the  scattering  of 
a periodic  acoustic  source  with  the  two-zone  overlap  configuration  shown 
schematically  in  Fig.  2a. 

The  single-domain  grid,  consisting  of  361  x 321  points,  is  split  along 
9 = 90°,  where  extra  £- lines  are  added  to  form  a five-point  overlap  as  in 
Fig.  1.  Solutions  are  advanced  separately  on  each  subdomain,  and  informa- 
tion is  exchanged  at  the  overlap  points  in  the  manner  previously  discussed 
in  Section  3.5.  The  C6  scheme  is  employed  for  interior  points  along  with 
fourth-  and  fifth-order  compact  operators  at  the  boundary  and  next-to- 
boundary  points  respectively  whereas  RK\  is  utilized  for  time-integration. 
In  the  interior  of  each  domain,  a lOth-order  filter  is  utilized  while  high-order 
one-sided  techniques,  described  in  Refs.  [3,  7],  are  invoked  near  boundaries. 
For  all  filter  operators,  the  coefficient  otf  = 0.45  is  specified. 

Figure  2b  displays  instantaneous  pressure  contours  in  the  vicinity  of 
the  cylinder.  It  is  apparent  that  the  pressure  waves  cross  the  grid  interface 
without  producing  any  noticeable  disruptions  of  the  interference  pattern 
even  though  pressure  waves  generated  by  the  source  propagate  through  the 
overlap  region  in  an  oblique  direction  to  the  zonal  interface.  A quantitative 
comparison  of  the  single-domain,  multiple-domain  and  analytic  solutions 
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is  given  in  Fig.  2c  in  terms  of  the  directivity  of  the  radiated  sound  at 
r/D  = 5.  The  directivity  obtained  with  a 6th-order  near-boundary  filter 
in  the  overlap  region  is  essentially  the  same  as  the  corresponding  single- 
domain baseline  solution,  and  both  results  are  in  excellent  agreement  with 
the  theoretical  solution.  These  results  highlight  the  potential  of  the  present 
high-order  methodology  for  domain-decomposition  applications  on  parallel 
computers.  Additional  details  on  the  manner  in  which  the  high-order  filter 
may  be  applied  to  yield  a robust  farfield  radiation  treatment  for  acoustic 
simulations  may  be  found  in  Ref.  [19]. 

4.2.  DNS  OF  A PULSED  WALLJET 

As  another  example,  the  method  has  been  employed  to  simulate  the  three- 
dimensional  transition  of  a forced,  finite  aspect-ratio,  plane  wall  jet  [9],  The 
wall  jet  configuration  considered  is  shown  in  Fig.  3a.  The  main  parameters 
governing  the  flow  are  the  Reynolds  number,  the  disturbance  characteristics 
at  the  jet  nozzle,  the  aspect  ratio  of  the  channel,  and  the  length  of  the  wall. 
In  the  present  study,  the  Reynolds  number,  based  on  jet  maximum  velocity 
{Umax)  and  nozzle  height  (h),  is  2150.  The  aspect  ratio  of  the  channel  is 
2b  jh  = 20.  The  mean  velocity  profile  in  the  normal  direction  at  the  nozzle 
exit  is  parabolic  and  corresponds  to  a fully-developed  laminar  channel  flow 
in  that  direction.  The  flow  is  forced  at  the  nozzle  exit  with  a frequency  of 
200 Hz  and  amplitude  corresponding  to  6%  of  the  jet  centerline  velocity. 

The  overall  flow  structure  is  shown  in  Fig.  3b  in  terms  of  an  iso-surface 
of  vorticity  magnitude.  The  transition  process  begins  with  the  formation 
of  shear-layer  and  wall  vortex  pairs  which,  due  to  the  energetic  forcing, 
appear  close  to  the  nozzle  exit  and  are  phased-locked  for  a short  distance 
downstream.  In  the  process  of  their  spanwise  evolution,  the  rollers  are  first 
split  into  a double-helical  structure,  which  is  clearly  discernable  near  the 
sidewalls.  This  feature  propagates  toward  the  center  while  also  expanding 
in  the  radial  direction.  The  spiral  vortex  branches  are  wound  in  a sense 
opposite  to  that  of  the  swirl  direction  of  the  vortex,  but  consistent  with 
the  direction  of  the  induced  axial  flow  which  exists  within  the  vortex  core. 
This  vortex  branching  and  helical  twisting  spreads  rapidly  through  self- 
induction  effects,  and  eventually  reaches  the  symmetry  plane  where  the 
vorticity  magnitude  within  the  vortex  core  is  drastically  diminished  (hence 
the  apparent  break  in  the  iso-surface,  region  ‘7’  in  Fig.  3b).  The  ability  of 
of  the  present  high-order  solver  to  discern  the  fine-scale  breakdown  to  tur- 
bulence is  shown  in  Fig.  3c  which  displays  contours  of  vorticity  magnitude 
on  a horizontal  plane  ( y/h  = 0.5).  Additional  details  of  the  simulation,  as 
well  as  comparison  with  high-resolution  experimental  measurements  may 
be  found  in  Ref.  [9]. 
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4.3.  LES  OF  A SPATIALLY-EVOLVING  SUPERSONIC  BOUNDARY  LAYER 

This  simulation  considers  a zero-pressure-gradient,  flat-plate  boundary  layer 
at  a Mach  number  of  2.25  and  Reg  ~ 6000,  which  corresponds  to  the  com- 
putational studies  of  Refs.  [20,  21]  Details  of  the  stretched  Cartesian  grids 
and  the  initial/boundary  conditions  employed  are  given  in  Ref.  [11].  In  this 
section,  only  results  (denoted  as  ‘no-model’)  obtained  using  the  lOth-order 
low-pass  filter  without  the  inclusion  of  an  SGS  model  are  considered.  So- 
lutions for  both  the  Smagorinsky  and  dynamic  subgrid-scale  stress  models 
can  be  found  in  Ref.  [11]. 

Figure  4a  shows  the  spanwised- averaged,  mean  skin- friction  coefficient. 
Downstream  of  the  transition  location  (which  is  sensitive  to  the  numeri- 
cal scheme  and  forcing  employed),  the  present  results  compare  favorably 
with  those  of  Ref.  [20].  This  fact  is  encouraging  since  a coarser  mesh 
(371  x 61  x 151)  is  used  in  the  present  computations  with  the  sixth-order 
compact  approach.  The  calculations  of  Ref.  [20]  were  performed  on  a much 
finer  grid  (971  x 55  x 321)  utilizing  a fifth-order  upwind-bias  algorithm  for 
the  convective  terms.  Reasonable  agreement  in  terms  of  the  mean  stream- 
wise  velocity  profile  is  also  shown  in  Fig.  4b.  Finally,  contours  of  the  instan- 
taneous spanwise  vorticity  component  at  a height  of  y+  ss  1.0  are  shown  in 
Fig.  4c  and  display  the  longitudinal  structures  characteristic  of  turbulent 
wall-bounded  flows. 

4.4.  BOUNDARY-LAYER  TRANSITION  OVER  A FLEXIBLE  PANEL 

As  a final  example  of  a simulation  of  multi-disciplinary  physics  with  the 
present  methodology,  consider  a transitional  boundary-layer  flow  over  a 
flexible  finite  panel  embedded  in  a rigid  surface  as  shown  schematically  in 
Fig.  5a.  This  problem  is  closely  related  to  classic  panel  flutter  phenomena, 
as  well  as  to  viscous  flow  over  compliant  surfaces.  The  panel  of  length  a 
and  thickness  h extends  over  the  region  0.5  < x/a  < 1.5.  The  leading-edge 
region  of  the  plate  (—0.5  < x/a  < 0.0)  is  formed  by  an  ellipse  of  half- 
thickness 0.05/i  (i.e.  aspect  ratio  10).  An  additional  challenge  posed  by  this 
aeroelastic  simulation  is  the  need  to  accomodate  the  surface  deflection  with 
a dynamically  deforming  mesh.  The  problem  has  been  examined  in  great 
detail  in  Ref.  [18]  which  should  be  consulted  for  details  regarding  boundary 
condition  implementation  and  mesh  resolution  studies. 

For  brevity,  only  select  results  obtained  at  M0 0 = 0.8  and  Rea  = 105 
are  summarized  here  to  highlight  the  ability  of  the  method  to  capture  the 
complicated  unsteady  phenomena  under  the  influence  of  flow-induced  sur- 
face deformation.  At  low  values  of  the  dynamic  pressure,  a steady  flow  is 
obtained  despite  the  adverse  pressure  gradient  induced  by  the  downward 
deflection  of  the  panel.  At  higher  dynamic  pressures,  however,  a travelling- 
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wave-flutter  phenomenon  is  observed  as  summarized  in  Fig.  5b,  c.  The 
instantaneous  panel  shapes  (not  shown)  display  a seventh-mode  oscillation 
with  a dominant  nondimensional  frequency  St  = fa/U <*,  = 1.52  which  is 
substantially  higher  than  the  fundamental  frequency  of  the  elastic  plate. 
The  high-mode  flexural  deflections  are  observed  to  travel  along  the  panel 
and  to  reflect  at  the  panel  edges.  These  high-frequency  fluctuations  re- 
sult in  a pronounced  acoustic  radiation  pattern  above  the  vibrating  plate, 
shown  in  Fig.  5b  in  terms  of  the  instantaneous  pressure  field.  Downstream 
of  the  flexible  surface,  a regular  train  of  vortical  disturbances  is  observed 
(Fig.  5c)  with  characteristic  wavelength  and  frequency  compatible  with 
those  of  Tollmien-Schlichting  (T-S)  instability.  The  travelling  wave  flutter 
appeared  to  originate  from  the  coupling  of  the  T-S  waves  with  the  panel 
high-mode  transverse  fluctuations,  and  this  convective  instability  ceases 
below  a critical  value  of  Reynolds  number. 
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Figure  1.  Schematic  of  mesh  overlap  for  multi-domain  strategy 
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Figure  2.  Multi-Domain  Acoustic  Scattering  Simulation 
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Figure  5.  Transitional  Boundary-Layer  Flow  Over  a Fluttering  Elastic  Panel 


